(library (grid-gen HD)

(export HD:vrs:rpuvw HD:vrs:rhuvw
        HD:vrs:q
        HD:vrs-proj
        HD:vrs-bkproj
        HD:q HD:F HD:chrs HD:L HD:R
        HD:qlen
        HD:vect
        HD:A
        HD:srcs
        HD:Potential
        HD:DS
)

(import (guile)
        (rnrs)
        (except (srfi srfi-1) map)
)
(import (grid-gen main))
(import (grid-gen vars))
(import (grid-gen phys))
(import (grid-gen math))
(import (grid-gen symath))
(import (grid-gen geom))

(define HD:qlen (lambda (nd) 5))

(define HD:vect (lambda (nd)
  (apply comp-vector
         (append (list (comp-float "rho")
                       (comp-float "p")
                 )
                 (list-head (list (comp-float "u")
                                  (comp-float "v")
                                  (comp-float "w")
                            )
                            nd
                 )
                 (list-head (list (comp-float "x")
                                  (comp-float "y")
                                  (comp-float "z")
                            )
                            nd
                 )
         )
  )
))

(define eps (lambda (vr)
     `(/ ,(vr 'p) (* (- ,PHYS:gam 1.0) ,(vr 'rho)))
))

(define Egas (lambda (vr)
  `(* ,(vr 'rho) (+ ,(vr 'eps) (* ,(vr 'u2) 0.5)))
))

(define ent->p (lambda (vr)
  `(* (/ (- ,PHYS:gam 1.0) ,PHYS:gam) 
      ,(vr 'rho)
      (- ,(vr 'h) (* ,(vr 'u2) 0.5))
   )
))

(define ent (lambda (vr)
  `(+ (* (/ ,(vr 'p) ,(vr 'rho)) (/ ,PHYS:gam (- ,PHYS:gam 1.0)))
      (* ,(vr 'u2) 0.5)
   )
))

(define HD:vrs:def (lambda (vl)
  (cond ((eq? vl 'rho) (int-err "Default rho access"))
        ((eq? vl 'p) (int-err "Default p access"))
        ((eq? vl 'u) 0.0)
        ((eq? vl 'v) 0.0)
        ((eq? vl 'w) 0.0)
        ((eq? vl 'h) (int-err "Default h access"))
        ((eq? vl 'E) (int-err "Default E access"))
        ((eq? vl 'eps) (int-err "Default eps access"))
        ((eq? vl 'c)   (int-err "Default c access"))
        ((eq? vl 'vl) (int-err "Default vl access"))
        ((eq? vl 'vr) (int-err "Default vr access"))
        (#t (format #t "C: ~a\n" vl)
            (int-err "Default ? access")
        )
  )
))

(define HD:vrs:rpuvw (lambda (rho p u v w)
  (letrec* ((self (lambda (vl)
    (cond ((eq? vl 'rho) rho)
          ((eq? vl 'p) p)
          ((eq? vl 'u) u)
          ((eq? vl 'v) v)
          ((eq? vl 'w) w)
          ((eq? vl 'u2) `(+ (** ,u 2.0) (** ,v 2.0) (** ,w 2.0)))
          ((eq? vl 'h)   (ent  self))
          ((eq? vl 'E)   (Egas self))
          ((eq? vl 'eps) (eps  self))
          ((eq? vl 'c)  `(sqrt ,(self 'c2)))
          ((eq? vl 'c2) `(/ (* ,PHYS:gam ,p) ,rho))
          ((eq? vl 'temp) `(/ ,(self 'p) (* ,(self 'rho) ,PHYS:Rgas)))
      ;    ((eq? vl 'vl) self)
      ;    ((eq? vl 'vr) self)
          (#t (HD:vrs:def vl))
    ))))
    self
  )
))

(define HD:vrs:rhuvw (lambda (rho h u v w)
  (letrec* ((self (lambda (vl)
    (cond ((eq? vl 'rho) rho)
          ((eq? vl 'p) (ent->p (HD:vrs:rhuvw rho h u v w)))
          ((eq? vl 'u) u)
          ((eq? vl 'v) v)
          ((eq? vl 'w) w)
          ((eq? vl 'u2) `(+ (** ,u 2.0) (** ,v 2.0) (** ,w 2.0)))
          ((eq? vl 'h) h)
          ((eq? vl 'E)   (Egas self))
          ((eq? vl 'eps) (eps  self))
          ((eq? vl 'c)  `(sqrt ,(self 'c2)))
          ((eq? vl 'temp) `(/ ,(self 'c2) (* ,PHYS:gam ,PHYS:Rgas)))
          ((eq? vl 'c2)
          `(* (- ,h 
                  (* (+ (** ,u 2.0) (** ,v 2.0) (** ,w 2.0))
                     0.5
                  )
              )
              (- ,PHYS:gam 1.0)
           )
          )
       ;   ((eq? vl 'vl) self)
       ;   ((eq? vl 'vr) self)
          (#t (HD:vrs:def vl))
    ))))
    self
  )
))

(define HD:vrs:q (lambda (qr)
  (let* ((ql (arr-elements qr))
  
         (q1 (list-ref ql 0))
         (q2 (list-ref ql 1))
         (q3 (list-ref ql 2))
         (q4 (list-ref ql 3))
         (q5 (list-ref ql 4))
         
         (rho q1)
         (u `(/ ,q2 ,q1))
         (v `(/ ,q3 ,q1))
         (w `(/ ,q4 ,q1))
         (p `(* (- ,PHYS:gam 1.0)
                (- ,q5
                   (/ (+  (** ,q2 2.0) (** ,q3 2.0) (** ,q4 2.0))
                      (* ,q1 2.0)
                   )
                )
             )
         )
        )
    (lambda (vl)
      (cond ((eq? vl 'q1) q1)
            ((eq? vl 'q2) q2)
            ((eq? vl 'q3) q3)
            ((eq? vl 'q4) q4)
            ((eq? vl 'q5) q5)
        (#t ((HD:vrs:rpuvw rho p u v w) vl))
      )
    )
  )
))

(define HD:vrs-proj (lambda (vr xi geom)
  (let* ((vl (math-cvector (vr 'u) (vr 'v) (vr 'w)))
         (vp (geom:append 
               (list-head (arr-elements ((geom 'proj) vl xi)) (length xi))
               0.0
             )
         )
        )
    (HD:vrs:rpuvw  (vr 'rho)
                   (vr 'p)
                   (list-ref vp 0)
                   (list-ref vp 1)
                   (list-ref vp 2)
    )
  )
))

(define HD:vrs-bkproj (lambda (vr xi geom)
  (let* ((vl (math-cvector (vr 'u) (vr 'v) (vr 'w)))
         (vp (geom:append 
               (list-head (arr-elements ((geom 'bkproj) vl xi)) (length xi))
               0.0
             )
         )
        )
    (HD:vrs:rpuvw  (vr 'rho)
                   (vr 'p)
                   (list-ref vp 0)
                   (list-ref vp 1)
                   (list-ref vp 2)
    )
  )
))

(define HD:q (lambda (vr)
  (math-cvector (vr 'rho) 
               `(* ,(vr 'rho) ,(vr 'u))
               `(* ,(vr 'rho) ,(vr 'v))
               `(* ,(vr 'rho) ,(vr 'w))
                (vr 'E)
  )
))

(define HD:F (lambda (d)
  (lambda (v0)
    (let* ((vr (HD:vrs:swap v0 d))
           (F0 (math-cvector `(* ,(vr 'rho) ,(vr 'u))
                             `(+ (* ,(vr 'rho) (** ,(vr 'u) 2.0)) ,(vr 'p))
                             `(* ,(vr 'rho) ,(vr 'u) ,(vr 'v))
                             `(* ,(vr 'rho) ,(vr 'u) ,(vr 'w))
                             `(* ,(vr 'rho) ,(vr 'u) ,(vr 'h))
               )
           )
           (F1 (swap-str F0 1 2))
           (F2 (swap-str F1 2 3))
          )
      (list-ref (list F0 F1 F2) d)
    )
  )
))

(define HD:A0 (lambda (vm)
  (math-array '(5 5)
              `(0.0
                1.0
                0.0
                0.0
                0.0
                
                (- (* (* 0.5 (- ,PHYS:gam 1.0)) ,(vm 'u2)) (** ,(vm 'u) 2.0))
                (* (- 3.0 ,PHYS:gam) ,(vm 'u))
                (* -1.0 (- ,PHYS:gam 1.0) ,(vm 'v))
                (* -1.0 (- ,PHYS:gam 1.0) ,(vm 'w))
                (- ,PHYS:gam 1.0)
               
                (* -1.0 ,(vm 'u) ,(vm 'v))
               ,(vm 'v)
               ,(vm 'u)
                0.0
                0.0
                
                (* -1.0 ,(vm 'u) ,(vm 'w))
               ,(vm 'w)
                0.0
               ,(vm 'u)
                0.0
                
                (* ,(vm 'u)
                    (- (* 0.5 (- ,PHYS:gam 1.0) ,(vm 'u2)) ,(vm 'h))
                )
                (- ,(vm 'h) (* (- ,PHYS:gam 1.0) (** ,(vm 'u) 2.0)))
                (* -1.0 (- ,PHYS:gam 1.0) ,(vm 'u) ,(vm 'v))
                (* -1.0 (- ,PHYS:gam 1.0) ,(vm 'u) ,(vm 'w))
                (* ,PHYS:gam ,(vm 'u))
               )
  )
))

(define HD:vrs:swap (lambda (v0 d)
  (let* ((v1 (HD:vrs:rpuvw  (v0 'rho)
                            (v0 'p)
                            (v0 'v)
                            (v0 'u)
                            (v0 'w)
             )
         )
         (v2 (HD:vrs:rpuvw  (v0 'rho)
                            (v0 'p)
                            (v0 'w)
                            (v0 'u)
                            (v0 'v)
             )
         )
        )
    (lambda (vl)
      (cond
        ((eq? vl 'vl) (HD:vrs:swap (v0 'vl) d))
        ((eq? vl 'vr) (HD:vrs:swap (v0 'vr) d))
        ((eq? vl 'c)  (v0 'c))
        ((eq? vl 'c2) (v0 'c2))
        (#t ((list-ref (list v0 v1 v2) d) vl))
      )
    )
  )
))

(define HD:A (lambda (d)
  (lambda (v0)
    (let* ((vm (HD:vrs:swap v0 d))
           (A0 (HD:A0 vm))
           (A1 (swap-str (swap-col A0 1 2) 1 2))
           (A2 (swap-str (swap-col A1 2 3) 2 3))
          )
      (list-ref (list A0 A1 A2) d)
    )
  )
))

(define HD:R0 (lambda (vm)
  (math-array '(5 5)
              `(1.0
                0.0
                2.0
                0.0
                1.0
     
                (- ,(vm 'u) ,(vm 'c))
                0.0
                (* 2.0 ,(vm 'u))
                0.0
                (+ ,(vm 'u) ,(vm 'c))
               
               ,(vm 'v)
                (* 2.0 ,(vm 'c))
                (* 2.0 ,(vm 'v))
                0.0
               ,(vm 'v)
               
               ,(vm 'w)
                0.0
                (* 2.0 ,(vm 'w))
                (* 2.0 ,(vm 'c))
               ,(vm 'w)
               
                (- ,(vm 'h) (* ,(vm 'u) ,(vm 'c)))
                (* 2.0 ,(vm 'v) ,(vm 'c))
               ,(vm 'u2)
                (* 2.0 ,(vm 'w) ,(vm 'c))
                (+ ,(vm 'h) (* ,(vm 'u) ,(vm 'c)))
               )
  )
))

(define HD:R (lambda (d)
  (lambda (v0)
    (let* ((vm (HD:vrs:swap v0 d))
           (R0 (HD:R0 vm))
           (R1 (swap-str R0 1 2))
           (R2 (swap-str R1 2 3))
          )
      (list-ref (list R0 R1 R2) d)
    )
  )
))

(define HD:L0 (lambda (vm) ;;; WRONG!
  (mul2 `(/ (- ,PHYS:gam 1.0) (* 2.0 (** ,(vm 'c) 2.0)))
         (math-array '(5 5)
                     `( (+ (* ,(vm 'u2) 0.5)
                           (/ (* ,(vm 'u) ,(vm 'c)) (- ,PHYS:gam 1.0))
                        )
                        (* -1.0
                           (+ ,(vm 'u)
                               (/ ,(vm 'c) (- ,PHYS:gam 1.0))
                           )
                        )
                        (* -1.0 ,(vm 'v))
                        (* -1.0 ,(vm 'w))
                        1.0
                       
                        (* -1.0 (/ (* ,(vm 'v) (** ,(vm 'c) 2.0))
                                   (- ,PHYS:gam 1.0)
                                )
                        )
                        0.0
                        (/ ,(vm 'c) (- ,PHYS:gam 1.0))
                        0.0
                        0.0

                        (* -1.0 (/ (* ,(vm 'w) (** ,(vm 'c) 2.0))
                                   (- ,PHYS:gam 1.0)
                                )
                        )
                        0.0
                        0.0
                        (/ (** ,(vm 'c) 2.0) (- ,PHYS:gam 1.0))
                        0.0
                       
                        (- (/ (** ,(vm 'c) 2.0) (- ,PHYS:gam 1.0))
                           (* 0.5 ,(vm 'u2))
                        )
                       ,(vm 'u)
                       ,(vm 'v)
                       ,(vm 'w)
                       -1.0

                        (- (* ,(vm 'u2) 0.5)
                           (/ (* ,(vm 'u) ,(vm 'c)) (- ,PHYS:gam 1.0))
                        )
                        (* -1.0
                           (- ,(vm 'u)
                               (/ ,(vm 'c) (- ,PHYS:gam 1.0))
                           )
                        )
                        (* -1.0 ,(vm 'v))
                        (* -1.0 ,(vm 'w))
                        1.0
                      )
         )
  )
))

(define HD:DS0 (lambda (vm)
  (let* ((vl (vm 'vl))
         (vr (vm 'vr))
         (c   (vm 'c))
         (c2  (vm 'c2))
         (rho (vm 'rho))
         (drho `(- ,(vr 'rho) ,(vl 'rho)))
         (dp   `(- ,(vr 'p)   ,(vl 'p)))
         (du   `(- ,(vr 'u)   ,(vl 'u)))
         (dv   `(- ,(vr 'v)   ,(vl 'v)))
         (dw   `(- ,(vr 'w)   ,(vl 'w)))
        )
    (math-cvector
     `(/ (- ,dp (* ,rho ,c ,du)) (* 2.0 ,c2))
     `(/ (* ,rho ,dv) (* 2.0 ,c))
     `(/ (- (* ,c2 ,drho) ,dp) (* 2.0 ,c2))
     `(/ (* ,rho ,dw) (* 2.0 ,c))
     `(/ (+ ,dp (* ,rho ,c ,du)) (* 2.0 ,c2))
    )
  )
))

(define HD:DS (lambda (d)
  (lambda (vm0)
    (let* ((vm (HD:vrs:swap vm0 d))
           (DS0 (HD:DS0 vm))
           (DS1 (swap-str DS0 1 2))
           (DS2 (swap-str DS1 2 3))
          )
      ;(list-ref (list DS0 DS1 DS2) d)
      DS0
    )
  )
))


(define HD:L (lambda (d)
  (lambda (v0)
    (let* ((vm (HD:vrs:swap v0 d))
           (L0 (HD:L0 vm))
           (L1 (swap-col L0 1 2))
           (L2 (swap-col L1 2 3))
          )
      (list-ref (list L0 L1 L2) d)
    )
  )
))

(define HD:chrs (lambda (d)
  (lambda (vm)
    (let* ((u (list-ref (list (vm 'u) (vm 'v) (vm 'w)) d)))
      (math-cvector `(- ,u ,(vm 'c))
                     u
                     u
                     u
                    `(+ ,u ,(vm 'c))
      )
    )
  )
))

(define HD:Potential (lambda (x y z)
  0.0
))

(define HD:srcs (lambda (vr xc tau)
  (let* ((vrs `(,(lamcon 'X) ,(lamcon 'Y) ,(lamcon 'Z)))
         (pf (apply PHYS:Potential vrs))
         (grd1 (map (lambda (v) (simp-ex (MATH:simplify (diff-eq pf v))))
                     vrs
               )
         )
         (vl (let ((h (make-hash-table)))
               (for-each
                 (lambda (x) (hash-set! h (car x) (cadr x)))
                 (zip vrs xc)
               )
               h
             )
         )
         (grd2 (map (lambda (x) (subst-vars x vl)) grd1))
         (vl2 (let ((h (make-hash-table)))
               (for-each
                 (lambda (x) (hash-set! h (car x) (cadr x)))
                 (zip vrs (list 'x 'y 'z))
               )
               h
             )
         )
         (dbg (map (lambda (x) (subst-vars x vl2)) grd1))
        )
  ;;  (format #t "\nDBG1: ~a\n" (infix-out (car dbg)))
  ;;  (format #t "\nDBG2: ~a\n" (infix-out (cadr dbg)))
  ;;  (format #t "\nDBG3: ~a\n" (infix-out (caddr dbg)))
    (HD:vrs:rpuvw (vr 'rho)
                  (vr 'p)
                 `(- ,(vr 'u) (* ,tau ,(car   grd2)))
                 `(- ,(vr 'v) (* ,tau ,(cadr  grd2)))
                 `(- ,(vr 'w) (* ,tau ,(caddr grd2)))
    )
  )
))


)

(import (grid-gen main))

(set! PHYS:vrs:rpuvw HD:vrs:rpuvw)
(set! PHYS:vrs:rhuvw HD:vrs:rhuvw)
(set! PHYS:vrs:q HD:vrs:q)
(set! PHYS:vrs-proj HD:vrs-proj)
(set! PHYS:vrs-bkproj HD:vrs-bkproj)
(set! PHYS:q HD:q)
(set! PHYS:F HD:F)
(set! PHYS:chrs HD:chrs)
(set! PHYS:A HD:A)
;;; (set! PHYS:L HD:L) ;;; WRONG!
(set! PHYS:R HD:R)
(set! PHYS:DS HD:DS)
(set! PHYS:qlen HD:qlen)
(set! PHYS:vect HD:vect)
(set! PHYS:srcs HD:srcs)
(set! PHYS:Potential HD:Potential)
